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Abstract 

We present a numerical method for solving Weyl's embedding prob- 
lem which consists of finding a global isometric embedding of a positively 
curved and positive-definite spherical 2-metric into the Euclidean three 
space. The method is based on a construction introduced by Weingarten 
and was used in Nirenberg's proof of Weyl's conjecture. The target em- 
bedding results as the endpoint of an embedding flow in R'^ beginning 
at the unit sphere's embedding. We employ spectral methods to handle 
functions on the surface and to solve various (non)-linear elliptic PDEs. 
Possible applications in 3 -I- 1 numerical relativity range from quasi-local 
mass and momentum measures to coarse-graining in inhomogeneous cos- 
mological models. 

1 Introduction 

It is a classic result in differential geometry that for every 2-surface of topol- 
ogy, equipped with a positive-definite metric whose curvature is positive, a global 
isometric embedding into the Euclidean three space can be found. Moreover, 
the embedding is unique up to the isometries of the target Euclidean space 
which are rigid motions: rotations and translations as well as reflections [T5]. 
This was conjectured by Hermann Weyl and later proven by Alexandrov and 
Pogorelov [2 [T7] , and independently by Nirenberg [IS] . 

Isometric embeddings into flat space have a wide range of applications in 
general relativity. For a given isometric embedding in a curved ambient space 
they provide a reference surface, thereby fixing the intrinsic geometries. From 
the difference between the extrinsic geometries it is possible to define measures of 
quasi- local mass like the Brown- York and Kijowski-Liu-Yau masses pTl [3l l4l [14] . 
For more examples see Szabados' review article |19| . Such quasi-local mass mea- 
sures - not yet applied in (3-1-1) numerical relativity - could potentially be useful 
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to determine mass and momenta in binary black hole simulations. In addition, 
the uniqueness of the reference shapes makes them ideal for visualising 2-metrics. 
Another application can be found in 3+1 numerical cosmology, where isometric 
cmbcddings of S'^ surfaces are the principle step in a method for coarse-graining 
(averaging) of expansion and shear in inhomogeneous cosmological models pro- 
posed by MK in [T1[T1]. 

The isometric embedding equations constitute a three-dimensional (3D) non- 
linear, coupled system of PDEs, for which standard numerical methods can not 
be applied and analytical solutions are impossible to find. This could be a 
reason why no practical application in numerical simulations exists. The prob- 
lem has already been addressed twice: by Nollert and Herold in [TB], and by 
Bondarescu, Alcubierre and Seidel [2]. In the algorithm presented in [16], the 
surface in question is triangulated and then a corresponding polyhedron in R'^ 
is constructed, whose edges have a length approximating those of the source 
triangulation. The method has a serious drawback: multiple polyhedra exist 
fulfilling the length constraints between neighboring points, most of them con- 
verge to unwanted non-regular (non-differentiable) surfaces. For the method in 
[2] the three unknown embedding functions are expanded into harmonic poly- 
nomials of maximal degree /^^^ and the embedding problem is turned into a 
numerical minimization problem of a 3(/„ax + l)'^ dimensional functional, which 
vanishes only if the embedding is isometric. This renders the method crucially 
dependent upon the minimisation algorithm used which, among other things, 
must steer away from the false minima of the functional. The minimization 
becomes increasingly difficult and computationally costly by increasing the di- 
mensionality of the parameter space. 

Our algorithm is based on the continuity method used in Nirenberg's proof 
and is not effected by the aforementioned problems. Its convergence to the right 
solution is guaranteed by Nirenberg's theorem. The target embedding results 
as the final surface of an embedding flow beginning with the unit sphere's em- 
bedding in R"^ . In each step of the flow the solution of the linearised embedding 
equation (LEE) allows one to continuously deform the surface from a round unit 
sphere to the final shape. This requires the conformal factor linking target and 
spherical round 2-metric, which we obtain as a steady-state solution of the Ricci 
flow. However we note that, unlike the previous two approaches, it is strictly 
limited to surfaces whose curvature is positive everywhere: Even if the isometric 
embedding exists for a surface with non-positive curvature, the algorithm can 
not be applied. 

We present a variety of numerical methods, among others, a pseudo-spectral 
parabolic evolution scheme to solve the various elliptic PDEs appearing in the 
problem that could be interesting for other purposes in numerical relativity. 

The paper is organized as follows: In the next section we provide the math- 
ematical background for the paper, as well as description of the linearized em- 
bedding flow we use. In the third section we discuss the technical details of 
the algorithm and its implementation; in the fourth, we present the results of a 
concrete numerical test case. We state the final conclusions in the fifth section 
and in the appendix, where we also quote relevant but more technical results 
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obtained in this paper. 



2 The mathematical background 

The Weyl problem can be described as foUows: Given a 2-surface C of spheri- 
cal topology, equipped with a positive-definite metric q (target metric), find a 
sufficiently regular embedding into the Euclidean space 

$ : C ^ R3 (1) 

which preserves the 2-metric g, i.e. 

q^<^*5 (2) 

where 5 is the standard metric on R'^ and $* denotes the standard puUback. In 
this paper we shall assume that the metric and the embedding arc smooth. 

It is known that the embedding always exists if the metric is sufficiently 
regular and its curvature is positive everywhere |15[ [9] . It is also known to be 
unique up to the rotations, translation and reflexions in the Euclidean space 
(see for example [TH] for a review of the rigidity results) . 

The algorithm to determine $, which we present in this paper, is based on 
the continuity method used in Nirenberg's proof of existence and introduced by 
Weyl. It consists of three steps which we will now describe briefly, leaving their 
detailed discussion until the next subsections. 

At first wc compute the conformal factor relating the target metric to the 
metric of the round sphere with radius one (round metric) which is a steady- 
state solution of the Ricci flow. 

\ = e-^^q (3) 

It is known that the Ricci flow uniformizes the metric [S] , i.e. the flow converges 
to a constant curvature metric q ^ % whereby we obtain a. In general, the 
round metric and a are then given in arbitrary coordinates. This has to be 
corrected through a transformation, for which the coordinate representation of 
% and its embedding $o into R'^ take a standard well-known form. 
Then we construct a one-parameter family of metrics 

%^^l{t,aY% (4) 

where ri(t, a) is a function chosen such that is the target metric and the 
round metric, which allows one to morph one metric smoothly into the other. 

Finally, we perform the embedding flow: Beginning at i = the known 
standard embedding $o of the round metric is deformed in 'small' steps such 
that the induced metric of the deformed surface at each step matches \ until 
the target metric is reached. This is accomplished by solving the linearized 
embedding equations, which relate a small change of the metric tensor with a 
small deformation of the embedding functions. 

We now proceed to describe the three steps in detail. 
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2.1 The Ricci flow 



In the first step of our method we have to determine the conformal factor be- 
tween the target and round sphere metric. In Riemannian 2-manifolds the Ricci 
flow reduces to a single equation for the conformal factor 

p = 7^b(r)]-(7^b(r)]), (5) 

where 7?,[p] is the Ricci scalar of the metric p and (72, [p]) its surface average wrt 
q included to keep the surface area of p bound. TVj}] is related to the Ricci 
scalar of the target metric in the following way 

n[p]=e'P{n[q]+2'i^p), (6) 

with 'A being the Laplacian wrt to q. 

The flow converges for large r to a function p — T a for which er"^" q is a 
round sphere metric (see [SUH], the latter also in [S]). By simple rcscaling of a 
we can ensure that the sphere has the area of 47r. Now it is possible to construct 
a family of metrics joining target and round sphere metric 

\=\^{t,af. (7) 

where ri(t, a) should be chosen such that 'TV\;q\ is always positive. This is guar- 
anteed for cr) = e^'^*. We have also tested a different function to drive the 
embedding flow cr) = t[eF — 1) + 1, in which the change between consecutive 
embeddings in the flow, as explained later, is linear in t as well. 



2.2 Round metric in standard coordinates 

The steady-state solution of the Ricci flow is a round metric \ ~ e~^°' q in 
arbitrary coordinates. To proceed, a transformation to standard coordinates a;' 
is required for which its isometric embedding X* into R'^ takes the well-known 
form, unique up to rotations and reflections 

=x/r= sin 6' cos (y9 
=y/r= sin sin (/3 
X^ =z/r= cose, (8) 

where r = \/ + y"^ -\- z'^ and (0, ip) are spherical coordinates corresponding to 
x*. Note that the three functions on the right hand side of ([8]) constitute a 
real, orthogonal basis for the spherical harmonics with / = 1. This means that 
these functions, denoted as n*, are eigenfunctions of the Laplace operator on 
the sphere with eigenvalue —2 

An* = -2n\ i = l,2,3. (9) 
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They are also orthogonal and normalized in the sense that 

Air 

rii dA = — Jy , (10) 

where dA is the area element associated with the round metric. It is easy to 
check that the opposite is also true: Any three orthogonal and properly nor- 
malized functions satisfying ([9]) are related to (|8]) by an 0(3) linear mapping 
and therefore constitute an isometric embedding of the unit sphere themselves. 
Equation ([9]) in turn can be solved numerically even in non-standard coordi- 
nates. 



2.3 The linearized embedding equation 

Consider an S'^ surface I?, endowed with a coordinate system 9^ and embedded 
in by a mapping described by three functions X^(0^). The induced metric 
has the form of 

QAB - X\AX\g5,,. (11) 

If we deform the embedding by adding a small 5X^{9^)^ the metric changes 
according to 

5qAB = 2<5X\^X^,j)% (12) 

up to the linear order. Given the metric change 5q one can ask for the compat- 
ible deformation vector. Finding SX'^ involves solving ([T2|) . which is called the 
linearized embedding equations (LEE) . Through a variable transformation this 
linear system of three PDEs can be turned into to a single elliptic equation of 
the second order for Weingarten's "Verschicbungsfunktion" and two ODEs, see 
also [9l [15] for derivations. Let Y"^ denote the deformation vector field we seek 
and dij the metric deformation and let be the outward-pointing null normal. 
We decompose into the normal and tangential part: 

Y' = Ts' + P 

(13) 

and introduce new variables ua and w 

UA = Sk Y^A 
w = -e^^DaIb. 

e^^ denotes the area form and Da is the covariant derivative on the surface, 
y can be reconstructed from ua and w via 

T^A = UA - Kab 
DaIb = -^w(AB + ]j^dAB + ^ Kab, (14) 
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where Kab denotes the extrinsic curvature. On a surface of positive scalar 
curvature Kab is positive definite0 and therefore has an inverse {K~^)^^ . The 
variable ua is then related to the first derivative of w by 

UA = ^cab (cn - , (15) 

CD = -DaQbd e 
Finally w itself has to satisfy an elliptic equation 

Cw = T. (16) 
The elliptic operator C is defined as 

Cw ^ Da {{K-^^w^b) +ICw, (17) 
where /C = K'^aj and the right hand side of the equation is 

r = DA{{K-Y''cB)^K^'AdBce^''. (18) 

Equation (|16p . together with ((TS|) and is equivalent to the original LEE. 

C is self-adjoint with the standard scalar product {f,g) = f* g dA[q{t)] 
and thus has only real eigenvalues. Moreover, for any convex surface it has 
a three-dimensional kernel spanned by the components of the normal vector 
s^{9"^), s^{9^) and s^{9^) [9]. In case of a round sphere it is possible to demon- 
strate that £ = A -I- 2. 

Since C has a non-trivial kernel, it is not invertible. Nevertheless, equation 
(fTB|) has solutions if its right hand side is orthogonal to the kernel: 

(r, s') ^ 0. (19) 

The solution is unique up to adding a combination of the functions s*, i = 1,2, 3. 
Geometrically this ambiguity corresponds to the possibility of adding a rigid 
rotation generator to solutions of ([T^. 

3 Technical and numerical details 

In this section we explain the numerical and technical details of our implemen- 
tation. Our approach consists of three main steps: the Ricci flow to relate 
target and round metric, the / = 1 eigenvalue problem of the round metric 
Laplacian in non-standard coordinates, and the embedding fiow from the round 
metric's embedding to the target embedding. For each computational step we 
need to conduct a variety of high accuracy numerical operations on spherical 
surfaces (numerical integration, interpolation, (anti)-differcntiation, coordinate 
inversion, solving elliptic PDEs). For this reason spectral methods in com- 
bination with particular coordinates are the best choice as we explain in the 
following. 

^This is the precise reason why our method is limited to positively curved metrics. 
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3.1 Coordinate basis, polynomial basis, grid setup 

Instead of covering spherical surfaces with 2D coordinate maps we consider 2- 
surfaces as being embedded in some (fictitious) ambient Riemanian 3-space (for 
example the Euclidian space), equipped with some 3D quasi-Cartesian coordi- 
nate system {i*} and represent all surface tensors using this exterior coordinate 
basis. Polar coordinates {9, (j)) on the surface are merely used to label the grid 
points {9i,(j>j), i = 1, - ■ ■ ,Ng, j = 1, • • • , 2Ng. We use a Gauss-Legendre grid 
structure, the canonical grid, on which a surface integral of polynomials of de- 
gree 21^^^ = 2{Ng — 1) can be represented exactly by a finite sum. 

Our approach requires various numerical operations on the surface: function 
evaluations at non-canonical points ("eval"), numerical integration, function 
inversion ( "inv" ) , differentiation and anti-differentiation ( "int" ) . For this reason 
we chose to represented the shape function h and other surface tensors as an 
expansion in harmonic polynomials. 

h = Y^'^lhf"' + OiL^^ + 1) = X] ^M'" + Ol^max + 1), (20) 

Im Im 

where F'™ is the standard othonormalised basis and the other basis is defined 
by := (n'A/^'™)', where = x'^/f = {sm9 cos 0, sin ^ sine/), cos 0) is the 
radial unit normal. This basis is orthogonal wrt distinct Z-eigen-spaces. Afi is 
a list of constant complex null vectors that span the 2/ -I- 1 harmonics in each 
eigenspace. The null vectors are chosen as in [10], then both basis are related 
by a discrete Fourier transform in each eigenspace. 



3.2 Differentiation on the surface 

The simple form of the basis (n^A/"/™)' is practical for evaluation of functions 
off the grid. In addition its differentiatior|jis straightforward, for example the 
derivative in 3d coordinates takes the form of 

d,h = d,n^ ^ *[/i]'"(7iJ AAj")'-i/AA^™. (21) 

By placing din^ — ?> den^ in front of the sum above, one obtains deh etc.; another 
practical feature. 

We prescribe 2-surfaces through level set functions G / shape functions /i, 
see fig. ((TJ (left), in some ambient manifold (E,7,;j), as is common in numerical 
simulations of the Einstein equation in (3-1-1) dimensions. If a shape function h 
is given, the computational steps to calculate curvature tensors on the surface 
are as follows 

^Alternatively, the differentiation of an expansion in y'™s can be performed by evaluating 
di, dijY^"^ from di, dij^''^ . It is practical to tabulated the di, dijY''"^ initially. 
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where the expressions in blue are only necessary for solving the LEE. 

The computation of curvature tensors on the surface requires '-DiXi of tensors 
Xi, which might have a normal component Xs 7^ O7 as for example dih, 

iD,Xj ^ CD^X^y ~ K,,xs, (23) 

which is not simply the tangential part of the ambient derivative CDiXjY ■ As 
a consequence, the Laplacian of a function ip on C is calculated by 

«A^A = - T'd.ij - IC ds^, (24) 

where = T*j, g-''^ is this contraction with the ChristofFcl symbols of 7^-. 

3.3 Parabolic flow relaxation method 

The elliptic PDEs of type C{u) — V{u) ~ appearing in our approach are solved 
by considering the associated parabolic flow equation [^1 

u^C{u)~V{u) ''^^ u^{C{u)^V{u))'^=° (25) 

where ii is the time derivative of u and {•}^'~° := • — (•) removes a function's 
surface average {1 = mode). For linear elliptic operators C with a non-positive 
spectrum, except I = 0, arbitrary initial data evolves to a steady-state solution 
(at ii, ii, • • • = 0) of the parabolic PDE which is automatically a solution to the 
elliptic equation. The Ricci flow is a parabolic PDE with a non-linear elliptic 
part. The existence of its steady states and convergence to them is known [6l[8]. 

The parabolic PDEs we deal with are solved by means of a pseudo-spectral 
scheme, where we use the method of lines, 2. or 4. order Runge-Kutta (RK2, 
RK4) and spectral finite difference methods, see eq. (f2T|) . until \u\ falls below 
machine/spectral precision. Since we are not interested in the details of the 
solution during the relaxation period, but in a quick fall-off, we pick a coarse 
time-step At and RK2 (Heun), where we take the CFL-condition of the heat 
equation with explicit Euler At < j/(Aa;)^, v = 0.25 as an orientation. 

The algebraic operations to compute the rhs produce an aliasing error. 
Therefore, we filter modes > ?„,,,^ form the rhs. 

3.4 Computational steps: Ricci flow, EVP, LEE 

The computational steps of the main component 'Ricci flow' of our implemen- 
tation arc the following. A shape function h{x) and an ambient geometry, see 
fig. ([1]) (left), are required as input 

^ In principle, the spectral decomposition of functions on the surface would allow us to turn 
non-linear / linear elliptic PDEs into algebraic / linear systems of equations to be solved with 
Newton-Raphson's / splitting (relaxation) methods which are numerically more efRcient but 
more stringent with regard to the initial guess and the conditioning of involved matrices. The 
parabolic flow relaxation method on the other hand covers non-linear / linear elliptic PDEs 
and is more robust. Moreover, efficiency is not a concern for the 2D PDEs we deal with. 
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Figure 1: Left: Level surface {C,qij) given by G(5:*) = f — h{x^) at G = 
and induced metric qij in the Riemanian manifold (I],7ij). Right: Embedding 
flow at t = ti,t2,t3 in the Euclidean space. Shift vector field at each step 
induces a new embedding F : Xj*-^^ i-> X^*2) = ^(i) + Y^ that is in general off the 
canonical grid (dots). 

h{x) cq. (122]) evolve: CT = {7e[e2'^g]}/'=° a{x),°q{x)^j, 

where n[e'^''q] = e^" {n[q] + 2 'JAcr) and «Acr is computed as in eq. The 
spherical round metric "g^^ = e~^°'qij with TZ[^q] = 2 is the steady-state solution 
of the Ricci flow. It is given in non-standard coordinates in which its Laplacian 
takes the form °A» ~ e^°'''A». As a standard coordinate basis, we take the 
three orthonormalised / = 1 cigcnfunctions of "A, sec eq. ^ 

a{x) — >■ evolve: = {"An* + 271*}/'"" (x) n^{x) cr(x) 

in which % takes the well known simple form = 5ij — niUj. Identifying 
x^ with Cartesian coordinates in maps "q.y , a to the Euclidean space where 
the embedding flow is performed, sec fig. ^ (right). We construct a family 
of metrics \ij ~ fl{t,a)^ ■ such that = and = q with IZ^q] > 
and divide the time interval into N steps, i.e. t = 0,ti,t2, ■ ■ ■ |f|w hereby we 
slowly deform *g into the target metric so that the difference 

:= q{tn+i\, - (26) 

is small, where is the induced metric wrt /i(„), sec fig. ([T]) (left), at the 

nth step. We initialize the fiow with the shape function /i(o) = 1 corresponding 
to the standard embedding of the unit sphere. At every step wc compute 

h{n} eq. (HI]) compute: '-"^dij -J> evolve: w = {{C{w) - t) / K.}* 
w Y'\p "(„+i)(2;|p) ?^'(2;(«+l)|p') iY\djY\a,q^j)\Q 

*The number of steps and the values tn are arbitrary and depend on how far q is off the 
round metric, usually 3 < N < 7. 
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Figure 2: (6', (/))-contour plots of the scalar curvature 'R-^q\ in relaxation flow at 
i = and t = 0.07. 



where * denotes projecting out the lowest eigenvalue of C (see Appendix). The 
complications arise from the fact that the new embedding induced by Y"^ is 
shifting the points off the canonical grid, moreover, the target metric i^q and a) 
has to be transported under the mapping X^^^^j^^ = ^(„) + from one surface 
to the other. 

After N equi-distant steps /i(7v) the difference (pS)) between the target metric 
and the induced metric wrt /i(7v) is of the order N ■ 0{d). But this solution can 
be refined by recursively applying [Jhe above computational steps with 
fixed in eq. (|26p whereby dij converges to zero exponentially. 



4 Numerical test case 

In this section we apply our procedure to a numerical test case, where we pre- 
scribe a test shape function in the Euclidean space given by the function below. 
We solve the embedding problem numerically and compare the original embed- 
ding with the resulting one. More specifically, we compute the induced 2-metric, 
perform the Ricci flow to the round sphere, solve the EV problem and perform 
the embedding flow from the round sphere to the target shape, solving the 
linearized embedding equation at each step. 

As a test case we pick a cigar-shaped function (for illustration purposes) and 
add a randomized non-polynomial part 

h = ci{\ + C2Re[{Uinj)f + ^1 + C3i?e[(A/'|nj)2 + [Mln^f]^ 

ci,2,3 = 0.458, 0.5, 0.15 (for which the area Awl and 7?. > 0, /C > 0) and 7Vi,2,3 
[flare three randomly oriented null vectors. This way all Ira modes are occupied 
in a somewhat random manner. The surface defined by this shape function has 

^If necessary this refinement process can be applied at any intermediate time tn, such that 
the embedding flow does not move to far off the conformal metric flow. 
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Figure 3: Left: Exponential fall-off of "\/ /i„(7^) for n = 2, 3, 4 at three different 
resolutions Ng ~ 16,20,24. Right: "old" (^, <?!))-coordinates given in the new 
{9, (/))-coordinates. 

the following area and central n-moments of the curvature " •\/ /in (72.): 

A/ At: = 0.999285, ^^Ai2(^) = 0.280044, 

3^^*3(^= 0.197685, t^iiTV) = 0.342832, 

where 

;.„(.) :=(((.)-.)"), 

and iZ ~ cTZ with (■^) = 1. As shown in fig. ^ (left) the scalar curvature 
significantly differs from that of a round sphere initially, where cr = is set, 
but smoothes out towards 7?, « 2 in the Ricci flow (right), where we employ the 
relaxation method as explained in the last section. The fall-off is exponential, 
see fig.© (left), until it reaches a plateau of truncation error which converges 
to zero by increasing the spherical resolution. The second order scheme is com- 
putationally more efficient, since the parabolic flow is stable up to = 0.11 for 
RK2-Heun and up to = 0.16 for RK4 [^(measured for Ne = 16, dx = 0.186). 

During the Ricci flow the surface area ( A/{4:Tt) = 0.999285 ) is bound but 
not fixed through ea. ipSj) . The resulting round sphere metric e~'^°'qij has the 
total area oi A/ (An) = 1.007889, which we normalise to unity. In the next step, 
we compute the I = 1 eigenfunctions of this metric, as explained in Sec. 
13.41 As initial data for the relaxation method we picked n-'|t=o = where 
the parabolic flow was stable up to = 0.23 for RK2. The three steady-state 
solutions are orthonormalised through the Gram-Schmidt process and inverted 
n^n) — W{n) using Newton- Raphson's method (since we have access to diU^ 
). In fig. ([3]) (right) the old spherical coordinate lines fi-' (n) are shown on the 
new grid on which we re-evaluate a / the spherical round metric ^qij- The 
round metric in the canonical coordinates can be embedded into Euclidean 

Re[J\fi,2.3] = l/V434{-9, 17, -8), 1/V3786(-1, -43,44), l/\/78(7, -2, -5) and 
/m[M,2,3] = -l/\/3(l,l,l). 

^ On the same surface solving the heat equation we got f = 0.20 RH2-Hcun, u = 0.28 
RH4. 
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Figure 4: Left: Contours of difference — between consecutive shape 

functions /i(„) for n = 1,2,3 in embedding flow (black 3-2, blue 2-1, light blue 
1-0). Right: Embedding flow at t„=o.2,3 = 0.0,0.66,1.0 (outer-transparent, 
inner-transparent, inner-solid). 

space by the shape function ho = 1. During the embedding flow the shape is 
stepwise 'deformed' into the target shape function by solving the LEE at each 
step. Again the relaxation method with RK2 for v = 0.2 is employed. 

We choose a conformal factor that is linear in t then the change of /i(„) 
during the first steps of the flow n = 1,2,3, tn=oa,2.3 = 0.0,0.33,0.66,1.0. It 
is about linear in t as well which can be seen in fig. Q (left); the difference 
between consecutive shape functions is approximately constant along a fixed 
direction. The increase/decrease is strongest in the directions where the target 
metric differs the most from the round metric. This is also apparent in fig. 
(jH) (left) and fig. ([6]), where the embeddings of q{t,a)ij at different states of 
the flow are shown. The drifting of grid points as illustrated in fig. ([IJ is not 
visualised but it can be seen how the unit sphere stepwise morphs into the target 
shape. The difference between /i(3) and the original shape function, see fig. ([6]), 
is about 10^^; in order to refine the solution further we recursively apply the 
LEE, whereby dij the difference between target metric and induced metric as 
well as the intrinsic curvature converges exponentially with increasing n and 
increasing resolution, limited by a plateau of truncation error, sec fig. ([7]). Note 
that due to the non-uniqueness of the embedding the final shape function is 
arbitrarily shifted and rotated wrt to the original one. In fig. ([7]) (left) the 
L2-norm in each Z-eigenspace is shown. It is independent of rotations but not 
of the shifts, which is clearly visible in the different I = 1 modes. 

5 Conclusion 

In this article we presented a new numerical algorithm to solve the Weyl prob- 
lem. The basic idea stems from the method of continuity which served as an 
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Figure 5: Embedding flow around unit sphere (transparent) at ti = 0.33 (left) 
and is = 1 (right) 
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flow. Right: Difference of the area A and UniJ^) between the induced and 
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approach in Nirenberg's proof of the conjecture and goes back to Weingarten 
and Weyl. The linearized embedding equations play a central role. They can 
be reduced to a single linear elliptic PDE by variable transformations. Solving 
the LEE stepwise for two 'nearby' metrics on a known embedding of one of the 
metrics allows one to 'slowly' deform an initial shape into the desired embedding 
(embedding flow). In doing so, it is necessary to link the target metric to an 
initial metric, whose embedding into is known. For this purpose the round 
metric is well suited which is a steady-state solution of the Ricci flow on 5^ that 
can be evolved from any metric. But it is then in general given in arbitrary 
coordinates. 

Hence, apart from solving the embedding equations, a number of additional 
obstacles appear, which complicates a numerical implementation of the method 
of continuity, namely we need to find a suitable coordinate system, solve various 
(non)-linear elliptic PDEs,(anti)-differentiate functions on the surface, handle 
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Figure 7: Left: L2-norm: y X]L=-i ^ original shape function h and 

hn=7 for Ng = 16. The difference for higher modes and Z = 1 is apparent. Right: 
Difference between the centered original and final shape function's fj,„ converges 
with increasing resolution and flow time down to trunction error. 



the drifting of grid points under the embedding flow or the 50(3, 1) freedom of 
the solution etc. Another challenge is to keep the total numerical error of all 
steps small. 

Through the use of ovcr-dctcrmincd quasi-Cartesian coordinates on the sur- 
face as well as spectral methods and the parabolic flow relaxation method we 
have shown how to overcome these technical difficulties altogether without loos- 
ing numerical accuracy. Note that the compatibility with the simulations in 
(3 -t- 1) numerical relativity is automatically given, since our approach requires 
as input an arbitrary shape function as well as a Riemannian 3-manifold, both 
of which generate an admissible 2-metric on the surface. This allows direct 
applications in binary black hole simulations to measure quasi-local mass or 
to investigate the large-scale dynamics of numerical simulations of inhomoge- 
neous cosmological models, where the coarse-graining process involves solving 
the isometric embedding problem. These quasi-local mass definitions, like the 
Brown- York or the Kijowski-Liu-Yau mass [TTl |31 IH [2] , are based on the com- 
parison principle: anchor the intrinsic geometries by isometric embeddings and 
compare the extrinsic geometries. Moreover, isometric embeddings are ideal for 
visualising 2-mctrics because they are fixed by unique shapes in the Euclidean 
space. 

We have implemented our method with Fortran90 and tested it on a concrete 
arbitrarily chosen initial shape function in the Euclidean space. By direct com- 
parison between the final and initial shape functions and metrics we were able 
to demonstrate exponential convergence by increasing the spherical resolution 
as well as the number of recursive iteration steps of the embedding flow. Our 
method is limited to positively curved and positive-definite 2-mctrics and to 
ones whose isometric embedding and shape function in R'^ is 'well represented' 
through harmonic polynomials within l^^^ < 40, i.e. whose spherical harmon- 
ics power spectrum drops quickly enough, but is then guaranteed to reach a 
solution at reasonable computational costs, a few minutes on a modern CPU. 

Apart from the applications already mentioned the numerical methods in the 
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sub-steps of our approach, like the implementation to solve the Ricci-flow or the 
I = 1 EV problem could lead to other applications in numerical relativity on their 
own. For example, to solve other elliptic PDEs that appear to find approximate 
Killing vector fields on trapping horizons [7] or to correct the gravitational waves 
signal of binary black hole simulations that is commonly extracted on coordinate 
spheres at finite radius. In general the spherical coordinates are distorted wrt 
the / = 1 eigenfunctions of these spheres. Furthermore, there are quasi-local 
mass and momentum measures, see |20] and references therein, which require 
isometric embeddings into Minkowski space (and into as a sub-step) that 
would allow one to measure the linear momentum in numerical BBH simulation 
quasi-locally. 

MK would like to thank Lars Andersson and Luciano Rezzola for their invi- 
tation to the Max Planck Institute for Gravitational Physics in Potsdam, which 
has enabled us to continue collaborating on this project. MJ thanks Badri Kr- 
ishnan, Alex Nielsen, Frank Ohme and Joachim Frieben for helpful comments 
and discussions. This work was supported by the IMPRS for Gravitational 
Wave Astronomy in the Max Planck Society. 

A Removing the positive eigenvalue from the 
spectrum of jC 

In the code we solve equation (fT6)) by relaxation method. This requires the 
linear operator C to be non-negative. This is not the case because of the positive 
second term in ((T7)) . C a positive principal eigenvalue Aq but can be modified 
£ to push that eigenvalue to the negative part of the spectrum. If the surface 
is not very far from a round sphere, we may expect that the spectrum of C will 
not be very different from the one of A -I- 2, and in particular will contain only 
one positive eigenvalue. In fact, our numerical evidence supports the conjecture 
that it is always the case for a convex surface. 

The eigenfunction corresponding to Ao is not known, so the modification of 
C is not completely trivial. First we define a new self-adjoint operator C as 

Vic \VicJ 

with /C denoting the trace of Kab (we assume Kab is negative definite). Since 
y/lC is positive everywhere, the dimensionality of subspaces where C is pos- 
itive, negative or zero is the same as in C. In particular, £ has only one 
positive eigenvalue, which is non-degenerate. One checks straightforward that 
C {VK^ — V^, so this eigenvalue and the corresponding function arc known 
explicitly. We now modify C further by defining 

Cig) = C{g) - , ^ . (Vic, g) VlC (29) 
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(28) 



for a constant C > 1. It has the only positive eigenvalue shifted to the negative 
part of the spectrum and thus can be tackled by the relaxation method. Namely, 
the parabolic equation 

/ = cif) - ^ + ^' (Vic, ^) Vic, (so) 



converges to the solution of ((T6)) divided by vX. We can simplify (pO)) further 
by introducing a new variable u = which then (j30[) becomes 

u=Y - ^) - / ^^^\ (1' ^(") - ^) ■ (31) 

^ (Vic, Vic) 



u tends exponentially to a solution of (|T6l). We have found out that this also 
holds true if we simplify this formula by replacing C / {ViC, ViC) by a sufficiently 
large, positive number. 

B Anti— differentiating a function in spherical har- 
monics 

In Section 12.31 we show that the vector field V can be reconstructed from its 
derivatives on the sphere. This leads us to the following numerical problem: 
Given the gradient of a function on expanded in terms of spherical harmonics, 
what is the spherical harmonics decomposition of the original function? 
Let / be a function expanded in terms of F'™ 

+00 m—l 

f = J2T. (32) 

;=o m=-l 

The (f) derivative is 

+00 m—l 

/.'^ = E E (33) 

;=0 m=-l 

SO for all coefficients with m 7^ we obtain a straightforward relation between 
the expansions of / and /^^i 

y[ff^ = -±y[fj^ . (34) 

m 

In general the 9 derivative of a regular / has a discontinuity at 6 ~ 0, n. This 
can be cured by multiplying it by sin6'. It is easy to see that sinfl^y'^™ can be 
expressed in terms of other spherical harmonics with the same m. In particular, 
the harmonics with m = satisfy 

sin4r'^o^^Ji±i^y'+- + ^J£±I^r'-^«. (35) 
oe ^(2Z + i)(2; + 3) V(2Z-i)(2; + i) 
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This equation gives us a recurrence relation for the m = coefficients of / in 
terms of sin6'^ 



y .o_ V(2?-l)(2^ + l) y W_1.0 
[/] Nn^/,,] 



+ /a + 1) V2[^ t""^^''^^ ^^^^ 

vahd for I > 2. For Z = 1 the relation reads 

n/]^'° = -^nsin0/,«]°^ (37) 

Finally the ^[/]*''° coefficient is arbitrary. This corresponds to the possibility of 
adding a constant to the solution. 



References 

[1] A. D. Alcxandrov. Existence of a polyhedron and of a convex surface with 
a given metric. Doklady Akademii Nauk S. S. S. R., 50:103-106, 1941. 

[2] Mihai Bondarcscu, Miguel Alcubierre, and Edward Scidel. Isometric em- 
bcddings of black hole horizons in three-dimensional flat space. Classical 
and Quantum Gravity, 19:375, 2002. 

[3] J.David Brown and Jr. York, James W. Quasilocal energy in general rel- 
ativity. 1991. In Mathematical Aspects of Classical Field Theory, edited 
by Mark J. Gotay, Jerrold E. Marsden, and Vincent Moncrief (American 
Mathematical Society, Providence, 1992) pages 129-142. 

[4] J.David Brown and Jr. York, James W. Quasilocal energy and conserved 
charges derived from the gravitational action. Phys.Rev., D47:1407-1419, 
1993. 

[5] H. D. Cao, B. Chow, S. C. Chu, and S. T. Yau, editors. Collected papers 
on Ricci flow, volume 37 of International Press Series in Ceometry and 
Topology. International Press, Sommerville, Massachusetts USA, 2003. 

[6] B. Chow. The Ricci flow on the 2-sphere. Journal of Differential Ceometry, 
33:325-334, 1991. 

[7] Gregory B. Cook and Bernard F. Whiting. Approximate Killing Vectors 
on S**2. Phys.Rev., D76:041501, 2007. 

[8] Richard Hamilton. The Ricci flow on surfaces. Mathematics and General 
Relativity, Contemporary Mathematics, 71:237-261, 1988. 

[9] Qing Han and Jia-Xing Hong. Isometric Embedding of Riemannian Mani- 
folds in Euclidean Spaces. American Mathematical Society, 2006. 



17 



[10] Michael Jasiulek. A New method to compute quasi-local spin and other 
invariants on marginally trapped surfaces. Class. Quant. Grav., 26:245008, 
2009. 

[11] Jerzy Kijowski. A simple derivation of canonical structure and quasi-local 
hamiltonians in general relativity. Gen. Relat. Grav. Journal, 29:307, 1997. 

[12] Mikolaj Korzyriski. Coarse-graining of inhomogcneous dust flow in General 
Relativity via isometric cmbcddings. AIP Conf. Proc, 1241:973-980, 2010. 

[13] Mikolaj Korzynski. Covariant coarse- graining of inhomogcneous dust flow 
in General Relativity. Glass. Quant. Grav., 27:105015, 2010. 

[14] Chiu-Chu Melissa Liu and Shing-Tung Yau. Positivity of Quasilocal Mass. 
Phys. Rev. Lett., 90:231102, 2003. 

[15] Louis Nirenberg. The Weyl and Minkowski problems in differential ge- 
ometry in the large. Gommunications on pure and applied mathematics, 
6:337-394, 1953. 

[16] Hans-Peter Nollert and Heinz Herold. Visualization in curved spacetimes. 
ii. visualization of surfaces via embeddings. In Relativity and Scientific 
Computing, pages 330-351. Springer- Verlag, Berlin, Heidelberg, 1996. 

[17] A. V. Pogorelov. On the proof of Weyl's theorem on the existence of a 
closed analytic convex surface realizing an analytic metric with positive 
curvature given on the sphere. Uspekhi Matematicheskikh Nauk, 32, 1949. 

[18] Michael Spivak. A Comprehensive Introduction to Differential Geometry, 
volume 5. Publish or Perish, Berkeley, 1979. 

[19] Laszlo B. Szabados. Quasi-local energy-momentum and angular momentum 
in gr: A review article. Living Reviews in Relativity, 7(4): 10, 2004. 

[20] M.-T. Wang and S.-T. Yau. Isometric Embeddings into the Minkowski 
Space and New Quasi-Local Mass. Communications in Mathematical 
Physics, 288:919-942, June 2009. 



18 



